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A  new  approach  is  developed  to  reduce  the  computational  complexity  of  a 
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dow,  a  traditional  batch  approach  would  result  in  a  large  number  of  multi¬ 
plication  and  add  operations  (i.e.,  an  order  N,  where  N  is  the  window  length). 
This  memorandum  shows  that  the  moving  average  batch  LMSF  procedure  could  be 
made  equivalent  to  a  recursive  process  with  identical  filter  memory  length  but 
at  an  order  of  reduction  in  computation  load.  The  increase  in  speed  due  to 
reduced  computation  could  make  the  moving  average  LMSF  procedure  competitive 
for  many  real-time  processing  applications. 
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1.  INTRODUCTION 

In  processing  a  long  stream  of  digital  data,  a  polynomial  moving  average 
Least  Mean  Square  Fit  (LMSF)  is  sometimes  used  to  provide  smoothing  or  filter¬ 
ing  of  the  noisy  component  imbedded  in  the  data  [1].  The  LMSF  has  a  finite 
duration  memory  which  Is  determined  by  the  length  of  the  moving  window. 

Figure  1  illustrates  the  successive  polynomial  LMSF  obtained  from  an  N-point 
window  which  slides  over  the  input  data  sequence.  These  operations  can  be 
further  explained  as  follows.  * 

Let  z(i)  for  i  =  1,  2,  ...  be  the  measurement  (data)  sequence,  and 
yn(k)  =  z(n  -  N  +  k)  for  k  =  1,  2,  ....  N  be  an  N-point  sequence  whose  first 
data  point  is  z(n  -  N  +  1)  and  the  last  data  point  is  z(n).  Note  that  we 

A 

assume  yn(k)  =  0  for  all  k  <  1.  Let  yn(k)  be  the  LMSF  estimate  based  on  the 
sequence  (y_(k);  k  =  1,  2,  ...,  N).  Then  for  a  given  k,  the  sequence  (y-i(k), 
YgU),  YgU),  ...),  obtained  from  a  repeated  LMSF  operations,  yields  a  smooth 
estimate  of  yn(k)  for  n  =  1,  2,  3,  ...  .  Note  that  for  1  <  k  <  N,  LMSF  per¬ 
forms  a  smoothing  operation,  and  for  k  =  N  the  LMSF  performs  a  filtering  func¬ 
tion  if  the  measurement  z(n)  is  made  available  at  the  current  time  [2].  Also 

a 

note  that  the  LMSF  estimates  yn(k)  and  yn+i(k)  contain  the  N-l  point  over¬ 
lapped  data  sequence  (z(n  -  N  +  2),  z(n  -  N  +  3),  ...,  z(n)).  Simply  stated, 
it  is  the  recognition  of  the  fact  that  moving  average  LMSF  computations  con¬ 
tain  the  overlapped  sequence  that  ultimately  allows  the  significant  saving  in 
computation.  The  reduction  in  computation  is  achieved  via  a  successful  devel- 

A  A 

opment  of  a  recursion  relation  between  yn(k)  and  yn+^(ic). 

There  are  several  advantages  in  using  a  moving  average  polynomial  LMSF 
[3].  Foremost  of  these  is  that  LMSF  provides  a  stable  operation  since  the 
LMSF  polynomial  coefficients  are  obtained  from  a  bank  of  Finite  Impulse 
Response  (FIR)  filters  (shown  later)  operating  on  the  data  sequence.  Sec¬ 
ondly,  the  finite  memory  window  assures  us  that  bad  data  points  which  lie  out¬ 
side  the  window  will  have  no  effect  on  the  resulting  LMSF  estimates. 

On  the  other  hand,  there  are  also  a  number  of  disadvantages  in  the  moving 
average  LMSF  application.  First  is  that  the  sizable  amount  of  memory  storage 
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641.172 


Figure  1.  Moving  Average  Least  Mean  Square  Fit 
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required  for  data  lies  within  the  operating  window.  Second  is  the  apparent 
large  computational  requirement  in  comparison  to  the  recursive  implementation 
using  Infinite  Impulse  Response  filters  [41.  The  study  presented  here  shows 
that  while  the  size  of  storage  requirement  is  unchanged,  the  computational 
load  can  be  reduced  significantly  via  an  equivalent  recursive  formulation. 

1.1  Problem  Statement 


Given  a  measurement  sequence  ( z ( i ) ) ;  i  =  0,  1,  2,  ....  we  denote  the 
first  N-point  sequence  by  the  column  vector 


Zy  =  (z(l),  z(2) . z(N))T 


(1) 


where  (  )  denotes  the  vector  transpose.  It  is  desired  to  find  a  vector  of 
Mth  order  polynomial  coefficients 


3|\j  -  (afj(0)j  <^(1),  •••» 


(2) 


that  minimizes  the  quadratic 


J(a*|)  "  -  H  -  H  a^,) 


(3) 


where  H  is  a  N  x  (M  +  1)  dimension  matrix  given  by 


H  = 


1  ll  l\ 
1  t2  t\ 


LM 


•  • 


•  • 


1  t  t^ 


tM 


(4) 


For  convenience,  we  shall  assume  a  uniform  sampling  interval  of  At 
seconds  such  that  one  can  write 


t.  =  (i  -  1 ) At  ;  i  -  1,  2,  3 


9  •  •  •  • 


(5) 
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Differentiating  Equation  (3)  with  respect  to  ^  and  setting  the  result 


to  zero  yields  the  desired  LMSF  coefficient  vector  given  by 


a.  =  (HT  H)"1  HT  z 


•=N 


=  (HT  H)"1  ^ 


(6a) 

(6b) 


where 


Xm  -  z.,  = 


-=-N 


1  1  1 


tx  t2  t3  ...  tN 


•  • 


,M  .  M  .M 
.  1  z2  z3  * 


lM 


■N  J 


N  N  N 

Ev,  •••  Z)*"2! 

,1=1  1=1  i=l 


(7a) 


(7b) 


Note  that  the  LMSF  polynomial 
M 

^(t)  =23  aN(m)tin  ; 


m=0 


0  <  t  <  NAt 


(8) 


can  be  used  to  interpolate  values  at  any  point  inside  the  data  window.  And  in 
particular  at  the  sampling  instants,  we  have 


^  =  H  ^  =  H(HT  H)'1  HT  z^ 


(9a) 


and  from  Equations  (1)  and  (9a)  we  obtain  for  the  kth  smooth  sample  (or  ele- 

a 

ment  in  z^) 


z(k)  =  (kth  row  of  H)^  ;  1  <  k  <  N 


(9b) 
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Also  note  the  subtle  differences  between  Equation  (6a)  and  Equa¬ 
tion  (6b).  First,  z^  is  an  N  dimensional  vector  whereas  is  a  M  +  1  dimen¬ 
sional  vector.  For  all  practical  purposes,  N,  the  number  of  data  points,  is 
much  larger  than  M,  the  order  of  the  LMSF.  Second,  the  matrix  (H^  H)~*  in 
Equation  (6a)  is  (M  +  1)  x  N  dimension,  and  the  matrix  (H^  H)"^-  in  Equa¬ 
tion  (6b)  is  (M  +  1)  x  (M  +  1)  dimension.  Thus  assuming  that  z^  and  x^  are 
known,  then  computing  a^  via  Equation  (6b)  will  result  in  saving  by  a  factor 
of  (M  +  1)/N  of  the  number  of  multiplications  and  adds.  For  example,  for  a 
typical  value  of  N  s  100  and  M  =  4,  we  have  (M  +  1)/N  =  .05.  This  is  a  sig¬ 
nificant  reduction  of  computation.  Of  course,  the  factor  (M  +  1)/N  is  the 
least  lower  bound  since  additional  computations  are  required  to  obtain  x^  from 
z^j.  If  Xj|  was  obtained  from  z^  directly  using  Equation  (7a),  then  no  improve¬ 
ment  is  gained.  Thus  it  is  desired  to  obtain  an  efficient  computation  of  x^ 
from  z^.  Consequently,  we  will  obtain  an  efficient  algorithm  to  calculate  ^ 
for  each  moving  average  window  of  length  N. 


Working  toward  this  goal,  we  generalize  Equation  (1)  in  writing 


*n 


(  (Z1’  z2’  zn>Ti  n±N 

(  (Zn-N+1’  zn-l*  zn^T  ;  n  > 


(10) 


as  an  N-point'  data  window  ending  at  the  sampling  instant  t  .  Then  the  LMSF 

A  I  I 

coefficient  vector  can  be  written  from  Equations  (6)  and  (10)  as 


In  =  (HT  H)-1  HT  ^ 


=  (HT  H)-1  Xh  ;  n  =  1,  2,  ...  . 


(11a) 

(lib) 


Note  that  given  the  smooth  sequence  ^  can  be  computed  from  the  relation 


4 =  H  a 


-n 


(11c) 


Therefore,  the  basic  statement  of  our  problem  is  to  obtain  an  efficient 

A 

methodology  to  calculate  a  for  all  n  via  obtaining  a  fast  recursion  in  com- 
puting  the  vector  2^.  The  matrix  (H  H)  is  specified  for  a  given  window 
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size  and  can  be  pre-computed .  Thus  matrix  inversion  is  not  required  for  each 
update . 

Note  that  letting  A  =  (H^  H)“*  in  Equation  (11a)  and  writing 


*n 


■  <sn<0>  V1' 


an(M)) 


Equation  (11a)  can  be  rewritten  as 


"yof 

■a„(l> 

• 

a 

• 

• 

A0(1 ) 
Ax(l) 


V1) 


Ao<2> 

Aj(2) 


W2> 


A0(N) 

'y„0)' 

• 

y„(2) 

• 

• 

• 

am(n). 

• 

• 

_y„(")_ 

(12) 


(13) 


Therefore,  the  (m  +  l)th  element  of  is 
N 

an(m)  =  2  V1)  yn^;  m  =  °» 1 . M  •  (14) 

i=l 

A  A 

From  Equation  (14)  we  note  that  each  element  of  a^,  say  an(m),  is  the 
output  of  a  Finite  Impulse  Response  (FIR)  filter  operating  on  the  data 
sequence.  Written  in  the  z-transform  notation,  the  FIR  filters  are  given  by 
N-l 

^  V1")  z1"N;  m  =  0,  1 . M  .  (15) 

i=0 


Equation  (15)  demonstrated  the  assertion  made  earlier  about  the  LMSF 
polynomial  coefficients. 

A 

A  final  point  we  want  to  make  is  that  the  LMSF  coefficient  an(m)  corre¬ 
sponds  to  the  mth  order  derivative  of  the  LMSF  polynomial  evaluated  at  time 
t  =  tx  =  0. 
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1.2  Technical  Objectives 


The  primary  objective  of  this  memorandum  is  to  present  the  step-by-step 
procedure  that  leads  to  the  derivation  of  the  fast-moving  average  least  mean 
square  recursion  algorithm.  The  performance  of  the  algorithm  will  be  verified 
through  extensive  computer  simulation.  Machine  timings  for  the  recursive 
formulation  will  be  obtained  and  used  to  compare  with  the  corresponding  batch 
approach . 

1.3  Report  Organization 

The  organization  of  this  report  is  as  follows.  Section  2  presents  the 
fast  moving  average  linear  least  square  recursion  algorithm.  Section  3 
extends  Section  2  results  to  the  general  least  square  case.  Section  4  dis¬ 
cusses  the  results  in  computer  simulations  and  timings.  Finally,  Section  5 
summarizes  the  results  presented  in  this  report. 


2.  MOVING  AVERAGE  RECURSIVE  LINEAR  LEAST  SQUARE 

In  this  section,  we  concentrate  on  the  recursive  algorithm  development 
to  the  linear  LMSF.  Partly  because  this  is  one  of  the  most  important  and 
widely  used  cases,  and  partly  because  the  study  on  the  linear  case  will  serve 
as  a  natural  base  for  extension  to  the  general  case. 

2.1  Problem  Formulation 


Let  (z { 1 ) ) ;  i  =  1,  2,  ...  be  the  measurement  sequence  and  t .  =  (i  -  l)At 
be  the  corresponding  sampling  time  where  At  is  the  sampling  interval.  A 
moving  average  first  order  fit  of  the  form 

z ( t )  =  aQ  +  aj  t  (16) 

with  an  N-point  data  window  is  desired.  We  seek  an  efficient  recursion  to 

A  A  *  T 

calculate  the  LMSF  coefficients  a^  =  (an(0),  an(l))  for  n  *  1,  2,  3,  ... 
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based  on  the  N-point  data  sequence  vector  ^  which  was  defined  earlier  in 
Equation  (10). 

2.2  Result  Derivation 

Suppose  the  first  N  data  points  are  available,  then  using  Equations  (4) 
and  (11) ,  we  obtain 

^  =  (HT  H)'1  HT  ^ 


(17) 


(18a) 


(18b) 
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N 

■  £ »,  *,  (iso 

i=l 

For  uniform  sampling  such  that  t.  =  (i  -  l)At,  the  following  quantities 
can  be  obtained. 

St,  -  11^^-  (19a) 

i=l 

.  at*  «(*  - 1) 

1=1  1 

D  .  At2  «2(sL^l) 
and  Equation  (17)  can  be  reduced  to 


(1%) 

(19c) 


yo> 

r  2(?n  -  }i 

6/At 

XNd) 

N(N  +  1) 

N(N  +  1) 

yn 

6/At 

12/At2 

xn(2) 

"  N(N  +  1) 

N(N2  1)_ 

(20) 


or  equivalently,  one  can  write 


B 


(21) 


where  the  matrix  B  -  (H^  H)"*  is  a  function  of  the  data  window  length  only.  We 
shall  next  develop  the  recursion  for  x„  as  follows. 


Let  n  =  N  +  l  where  l  is  an  Integer  such  that  for  any  given  integer  N  one 
can  write  n  =  1,  2,  ...  .  Then  from  Equation  (18b)  we  obtain 
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N 

XN+s/1)  =  ^2  Zi+Jt  *  (22a) 

i=l 

From  which  we  can  write 

N 

XN+JL^)  ‘  XN+A-1^^  *^2  (Z1+S.  ‘  Zi+S,-1^ 

i=l 

A 

=  ZN+fi,  "  zl  •  (22b) 

Substituting  n  =  N  +  l  in  Equation  (22b)  and  rearranging  terms  yields 
the  desired  recursion 

*p(l)  =  x^_^(l)  ^zn  ~  Zn-N^’  ^  -  1,  2,  ...  (22c) 

where  xQ(l)  -  0  and  zp_N  =  0  for  n  <  N. 

Now  from  Equation  (18c)  we  obtain 
N 

xN+fi/2)  =  ^  ti  Zi+S,  ^23a^ 

i*l 

and  so  for  t^  =  0,  we  obtain 

N 

xN+j,(2)  '  XN+A-1^2^  =  ^MZi+Jl  "  zi+£-l) 

i=l 

N-l 

^2  ^1  '  ti+l^Zi+S.  +  tN  ZN+S. 

1=1 

N 

"At  Xw  Zi+Jl  +  NAt  ZN+Jl 
1=1 

=  -At  Xj^_j_^(1)  +  NAt  z^^  .  (23b) 
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Again  substituting  n  =  N  +  JL  In  Equation  (23b)  yields  the  desired  recur 

sion. 


xn(2)  =  xnl(2)  +  NAt(zn  -  xn(l)/N)  (23c) 

where  xn(l)  is  obtained  from  Equation  (22c)  and  Xq(2)  -  0. 

2.3  Results  Analysis 

Thus  Equations  (22c)  and  (23c)  define  the  recursions  for  the  moving 
average  linear  least  square.  To  calculate  from  x^,  the  batch  approach 
requires  a  total  of  2N  sums  and  N  multiplications  (Equations  (22a)  and  (23a)), 
whereas  the  recursive  formulation  requires  only  four  sums  and  one  multiplica¬ 
tion  (Equations  (22c)  and  (23c)).  Note  that  the  latter  is  independent  of  the 
window  length. 

In  the  following  section,  we  will  extend  the  recursive  approach  to  the 
more  general  case. 


3.  MOVING  AVERAGE  RECURSIVE  GENERAL  LEAST  SQUARE 

In  this  section  we  generalize  the  results  obtained  for  the  linear  least 
square  problem  to  the  general  least  square  case.  Analytic  expression  will  be 
given  for  the  reduction  in  computation. 

3.1  Problem  Formulation 


Let  (z(i));  i  =  1,  2,  ...  be  the  measurement  sequence  at  time 
t^  =  (i  -  l)At .  A  moving  average,  N-point  window,  general  LMSF  using  poly¬ 
nomial  of  the  form 

z(t)  =  aQ  +  axt  +  a2t2  +  ...  +  aMtM  (24) 
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is  desired.  We  seek  an  efficient  recursion  to  calculate  the  LMSF  coefficients 

AAA  A  T 

in  =  ( an ( 0 ) ,  an(l)»  ...»  ap(M))  for  n  =  1,  2,  3,  ...  based  on  the  N-point 
data  sequence  vector  ^  defined  in  Equation  (10). 

3.2  Result  Derivation 


Suppose  the  first  N  data  points  are  available,  then  using  Equations  (4) 
and  (11),  we  obtain 

a,,  =  (HT  H)"1  HT  ^ 

=6)^  (25) 

where  B  -  (H^  H)“^  is  a  (M  +  1)  x  (M  +  1)  dimension  matrix  independent  of  the 
data  and 


■*N  =  hT  h 


1  1  1  ...  1 
tl  t2  t3  *  *•  *N 


•  •  • 


•  • 


•  •  t 


+M  +M 

x2  l3  * " 


N 

s*, 

i=l 

N 

St, 

i=l 


St;,, 

i=l 


(26) 
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Thus  the  mth  component  of  the  vector  is 
N 

V"1)  =  S  ti'1  zi  ;  m  -  1.  2 . M  +  1  .  (27) 

i=l 

We  want  to  develop  a  recursion  for  x^  for  n  =  1,  2,'  ...  . 

For  m  =  1  and  2,  the  recursions  are  identical  with  the  linear  least 
square  expressions  given  in  Equations  (22c)  and  (23c)  and  for  convenience 
restated  below  as: 

MU  -  *„-i<U  +  <*„  -  z„.n>  <28> 

xn(2)  =  x^jfZ)  +  N4t(zn  -  xn(l)/N)  (29) 

where  Xq(1)  =  Xq(2)  =  0  and  zn_^  =  0  for  n  <  N. 

Using  the  same  approach  as  in  the  linear  case,  the  recursion  for  the  gen¬ 
eral  case  (for  m  >'  2)  can  be  obtained  as  follows: 

N  N 

W«»  -  Vi-i<m)  ■  E**'1  zm  -  E^'1  zm-i 

1=1  1=1 


zi+fl,-l^ 


N-l 

iLi  (*1  1  '  t'l+l^Zi+il  +  zN+i, 
i=l 


um-l> 


-  t 


1  zs. 


(30) 


Since  t^  =  0,  tN  =  (N  -  l)At  and  using  the  binomial  expansion 
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(1  +  x)n  =  1  +  nx  +  ^  x2  + 

n 


...  +  x 


£  C<">  xJ 
j=o 


where 


(31a) 


r(n)  _  n! 

i  "  In  - 1 )  1  j  i 


(31b) 


is  the  binomial  coefficient,  we  can  write 


tf’1  -  tf;J  =  Atm-1  ((i  -  l)ra_1  -  i"1-1) 


=  Atn'“1  [(i  -  l)”*-1  -  (1  +  (i  -  l))1"-1] 


m-1 


At 


=  -At1 


m- 


1  (a  -  D™-1  -  2  ct1)  (1  - 1)J 
Vi  ja0 

m-i  ^  c(m-l)  (1  _  x)j  _ 

£0 


') 


(31c) 


Substituting  Equation  (31c)  in  Equation  (30)  and  simplifying  as  follows 
yields 

W")  -  *N+»-l<">  =  ££  cj""1^  -  1)J  Z1n  +  tj-1  ZN+J1 

i=l  j=0 


=  -At1 


m 


m-2  /  N-l  \ 

j=0  \  01  / 


-N+A 


m-2 

=  -At"1’1  £  Cj-U  At 
j=0 


Is *! 


7  .  7 

Zi+Jl  lN  N+i, 


)  +  trn'1 


N+A 


m-2 


=  _Atm-l  £  C(m-DAt-j  (xNH(j+l)  -  tj  zN+A)  +  tJ-1  zN+, 

j=0  v  ' 
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m-2 

■  %  w- E6!""11  wi*1)  (32a> 

j=o 

where  and  are  defined  by  the  relations 


m-1 


-  Dj  cf”-1' 

(32b) 

j*0 

e(m-l)  =  Atm-j-l  c(m-l)  > 

J  J 

(32c) 

Finally  letting  n  =  N  +  2.  in  Equation  (32a)  yields 

(m-2 
zn  ■  ^  S  8 

as  the  general  recursion  for  m  >  2  with  XQ(m)  =  0  for  all  m. 
3.3  Result  Analysis 


(m-1) 


(33) 


Equations  (28),  (29)  and  (33)  describe  the  recursion  for  computing^ 
given  x^.  Note  that  in  Equation  (33)  the  coefficient  oi^  and  can  be 

pre-computed  and  tabulated.  Note  also  that  in  order  to  compute  xn(m),  one 
must  first  compute  xn(l),  xn(2),  ...»  xn(m-l)  sequentially. 

The  following  development  briefly  examines  the  total  multiplications  and 
adds  required  for  the  recursive  approach.  A  study  on  the  computation  require¬ 
ment  in  Equations  (28),  (29)  and  (33)  shows  that  for  each  xn(m),  the  total 
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multiplications  and  adds  is  bounded  below  by  m.  Therefore,  for  an  Mth  order 
LMSF,  each  x_  update  requires  at  least  (M+l)(M+2)/2  multiplications  and  adds. 
To  update  the  polynomial  coefficient  vector  a^  there  are  additional  (M+l)^ 
multiplications  and  adds.  Thus  the  total  minimum  number  of  multiplications 

A 

and  adds  required  for  each  ^  update  is 

Nr  =  (M  +  l)2  +  (M  +  1) (M  +  2)/2 
=  (M  +  1) (M  +  1  +  (M  +  2)/2) 

=  (M  +  1 ) ( ( 3M  +  4)/2)  .  (34) 


However,  for  the  batch  LMSF,  it  was  shown  in  Section  1.1  that  the  correspond¬ 
ing  number  is 


Nb  =  (M  +  1)N 


(35) 


Therefore  the  recursive  approach  reduces  the  computation  load  with  respect  to 
the  batch  approach  by  a  factor 


Nr 


3M  +  4 

— zr~ 


(36) 


This  number  is  slightly  greater  than  the  ideal  lower  bound  (M  +  1)/N 
shown  in  Section  1.1 


4.  NUMERICAL  EVALUATION 

In  this  section,  we  discuss  the  method  of  implementation,  numerical 
accuracy,  computer  simulation  procedure,  and  the  result  of  computing  timings 
between  the  batch  and  the  recursive  approaches. 

4.1  Computational  Method 

Appendix  A  lists  the  computer  program  used  to  calculate  the  computa¬ 
tional  speed  between  the  batch  and  the  recursive  approaches.  The  following 
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discussions  briefly  highlight  some  of  the  implementation  consideration  and 
the  method  of  operation. 


To  prepare  for  the  recursion,  several  arrays  must  be  initialized  to  con¬ 
tain  the  constants  needed  in  the  calculations.  First  the  matrix  (H^  H)“*  is 
needed  where  the  H  matrix  is  the  matrix  previously  defined.  For  convenience 
we  assume  a  unity  sampling  interval  and  write  the  measurement  matrix  H(i,j)  = 
(i  -  1)J'*  where  the  row  variable  i  has  constraints  1  <  i  <  N  where  N  is  the 
length  of  the  data  window.  The  column  variable  j  has  constraints: 

1  <  j  <  M+l  where  M  is  the  fit  order  desired.  The  multiplication  of  HT  and  H 
is  trivial  but  taking  the  Inverse  creates  quite  a  numerical  problem  since  h"*"  H 
is  progressively  ill  conditioned.  A  single  precision  inversion  method  will 
only  provide  satisfactory  results  for  a  linear  or  first  order  fit;  double  pre¬ 
cision  will  yield  acceptable  results  for  second  and  third  order  fits,  but 
after  third  order  very  complicated  methods  will  be  necessary.  It  should  be 
pointed  out,  however,  that  the  matrix  inversion  needs  to  be  calculated  only 
once.  Thus  accurate  numerical  techniques  can  be  used. 


The  cjm) 


coefficients  are  now  calculated  using  the  formula 


m! 


(m  -  J 


where  M  is  the  fit  order  and  j  has  constraints  0  <  j  <  m  -  2.  Only  the  last 
M  -  1  rows  of  the  (M  +  1)  x  (M  +  1)  matrix  are  calculated  since  the  first  two 
rows  are  never  used  (see  Appendix  A). 

The  vector  is  calculated  next  using  the  formula  (see  Equation  (32b)) 


m-1 


%  ■  2Z (N  ■ 1)3 

jo 


r(m-l) 

Cj 


where  N  is  the  length  of  the  data  window  used. 
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Once  these  constants  have  been  prepared,  the  processing  loop  is  entered 
where  the  recursive  polynomial  coefficient  calculations  are  made. 

Data  is  first  read  in  and  then  the  y^  vector  (this  is  the  vector  x^  in 
Equation  (33))  is  calculated  using  the  cjm^  matrix  and  the  vector.  yN  is 
then  multiplied  by  the  (H^  H)"^  matrix  to  obtain  the  final  polynomial  coeffi¬ 
cients.  Next  an  estimate  of  the  next  path  point  is  calculated  using  the  first 
position  outside  of  the  data  window  as  the  target  point.  Once  this  is  done, 
the  oldest  data  point  is  discarded  and  is  replaced  by  the  newest  piece  of 
data.  A  pointer  moves  through  a  fixed  vector  to  keep  track  of  the  data  at  the 
back  of  the  data  window  (oldest  data  point)  so  that  it  can  be  used  in  the  next 
recursion.  The  last  step  is  to  make  the  current  recursion  results  "old"  by 
setting  the  yg (id)  vector  eclual  t0  the  ^^(ew)  vector*  This  loop  continues 
indefinitely  or  until  a  desired  condition  is  met. 


4.2  Numerical  Simulation 


The  testing  of  this  algorithm  was  done  on  DEC'S  VAX  11/780  located  at  the 
NUSC  Code  321  Advanced  Sonar  Development  Facility  (ASDF).  Fortran  IV  was 
used  in  its  programming.  For  uniform  testing,  the  same  set  of  10,000  data 
points  were  used  for  the  recursive  and  batch-moving  average  LMSF.  This  input 
data  consisted  of  a  unit  sampling,  linearly  increasing  ramp  with  unity  slope 
as  the  signal  and  an  additive,  zero  mean,  unit  variance  noise  sequence. 

Using  a  window  length  of  50  seconds  (i.e.,  N  =  50),  the  fitting  perform¬ 
ance  of  the  moving  average  recursive  LMSF  and  the  absolute  error  characteris¬ 
tics  as  a  function  of  time  are  presented  in  Figures  2  to  5  for  the  first, 
second,  third,  and  fourth  order  LMSF,  respectively.  There  are  two  important 
aspects  which  need  to  be  pointed  out.  First,  although  not  shown,  extensive 
testings  show  that  when  the  data  window  is  full,  both  batch  and  recursive 
approaches  yield  identical  results.  We  should  caution  the  reader  that  this 
may  not  be  the  case  when  a  numerical  problem  exists  in  matrix  inversion.  For 
example,  when  the  matrix  HT  H  Is  singular,  the  recursion  approach  fails  while 
there  exists  other  techniques  such  as  the  singular  value  decomposition  (SVD) 
in  solving  the  batch  approach.  Second,  prior  to  filling  the  data  window 
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TIME < SEC. > 

2a.  Recursive  Moving  Average  First  Order  LMSF  (Window  Length  =  50) 


Figure  2b.  Recursive  Moving  Average  First  Order  LMSF  Absolute  Error 

19 


ERROR <ABS  >  n  INPUT  DATA 


TM  841143 


TIME < SEC. > 

3a.  Recursive  Moving  Average  Second  Order  LMSF  (Window  Length  =  50) 


T I ME ( SEC . > 

Figure  3b.  Recursive  Moving  Average  Second  Order  LMSF  Absolute  Error 
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Figure  4a.  Recursive  Moving  Average  Third  Order  LMSF  (Window  Length  =  50) 


Figure  4b.  Recursive  Moving  Average  Third  Order  LMSF  Absolute  Error 
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TIME < SEC. > 

5a.  Recursive  Moving  Average  Fourth  Order  LMSF  (Window  Length  =  50) 


TIME  <  SEC .  > 

Figure  5b.  Recursive  Moving  Average  Fourth  Order  LMSF  Absolute  Error 
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(50  seconds),  the  recursive  approach  has  been  properly  Initialized  and  main¬ 
tains  relatively  good  fits  during  this  interval.  It  is  observed  that  higher 
order  LMSF  yield  smaller  fitting  errors  prior  to  the  data  window  being  fully 
filled.  On  the  other  hand,  the  first  order  LMSF  gives  the  best  result  after 
the  data  window  is  filled.  This  Is  to  be  expected  since  the  first  order  LMSF 
matches  the  input  signal  component.  Although  not  shown,  the  batch  approach, 
on  the  other  hand,  produces  no  output  at  all  before  the  data  window  is  full. 
Thus,  this  is  the  additional  advantage  for  the  recursive  approach. 

We  next  examine  the  relative  computation  time  required  for  the  batch  and 
the  recursive  approaches.  The  CPU  time  has  been  selected  as  the  performance 
measure  variable.  Using  a  number  of  different  length  windows,  both  approaches 
were  used  to  run  through  the  10,000  data  points  with  the  resulting  CPU  times 
properly  recorded.  The  results  are  presented  in  Figures  6  to  9.  Close  study 
of  these  results  shows  that  (1)  while  the  batch  CPU  time  increases  linearly 
with  the  window  length,  the  recursive  CPU  time  remains  constant  and  insensi¬ 
tive  to  the  window  length  as  expected,  (2)  for  a  given  window  length,  the 
batch  CPU  time  increases  (actually  linearly)  with  the  order  of  the  LMSF,  and 
(3)  as  window  length  shortens,  CPU  times  between  batch  and  recursive  becomes 
comparable. 

•4.3  Speed  Comparisons 

We  next  obtain  the  speed  ratio  (SR)  of  recursive  CPU  time  over  batch  CPU 
time.  The  results  are  shown  in  Figures  10  to  13  for  first,  second,  third,  and 
fourth  order  respectively.  Also  included  are  the  lower  bounds  obtained  based 
on  the  pure  number  of  multiplications  and  adds  required.  For  each  update  it 
was  shown  earlier  that  the  batch  approach  requires  approximately  (M  +  1)N 
number  of  operations.  For  the  recursive,  it  was  obtained  exactly  that  there 
are 


N.  =  l  +  (M  +  1)(3M  +  2) 
a  2, 


(37) 
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WINDOW  LENGTH < SEC . > 

Figure  6.  CPU  Time  Versus  Window  Length  Between  Batch  and  Recursive  First  Order  LMSF 


Figure  7.  CPU  Time  Versus  Window  Length  Between  Batch  and  Recursive  Second  Order  LMSF 
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0.  30.  100.  150.  200. 


WINDOW  LENGTH < SEC. > 

Figure  8.  CPU  Time  Versus  Window  Length  Between  Batch  and  Recursive  Third  Order  LMSF 


WINDOW  LENGTH < SEC. > 


Figure  9.  CPU  Time  Versus  Window  Length  Between  Batch  and  Recursive  Fourth  Order  LMSF 
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0.  50,  100.  150.  200. 


WINDOW  <  SEC .  > 

Recursive  Over  Batch  Speed  Ratio  Versus  Window  Length  (First  Order  LMSF) 


WINDOW < SEC . > 


Figure  11.  Recursive  Over  Batch  Speed  Ratio  Versus  Window  Length  (Second  Order  LMSF) 
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Figure  12.  Recursive  Over  Batch  Speed  Ratio  Versus  Window  Length  (Third  Order  LMSF) 


Figure  13.  Recursive  Over  Batch  Speed  Ratio  Versus  Window  Length  (Fourth  Order  LMSF) 
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number  of  adds  (Note,  M  is  the  order  of  the  fit)  and 


=  (M  ±  1)(3H  +  4) 
m 


1 


(38) 


number  of  multiplications.  Let  Ta  and  Tm  be  the  time  required  for  a  single 
add  and  multiplication,  respectively.  Then  the  lower  bound  (LB)  is  computed  from 
the  formula 


LB 


Na  Ta 

TW 


+  N 

m 

urr 


(39) 


Examining  Figures  10  to  13  shows  that  the  actual  SRs  fall  off  as  a  func¬ 
tion  of  window  length  at  a  rate  similar  to  the  lower  bound.  The  significant 
difference  between  the  actual  SR  and  the  LB  is  due  to  programming  overhead 
cost;  i.e.,  CPU  time  for  non-arithmetic  operations.  With  more  efficient  pro¬ 
gramming,  this  difference  is  expected  to  reduce.  At  any  rate.  Figures  10  to 
13  show  that  at  a  window  length  of  200,  the  recursive  CPU  time  is  only  20  per¬ 
cent  of  the  batch  CPU  time.  The  corresponding  number  for  the  lower  bound  is 
4  percent.  Note  also  that  in  contrast  to  the  lower  bound  behavior.  Figures  10 
to  13  show  that  for  a  given  window  length  N  the  actual  SRs  fall  off  as  a  func¬ 
tion  of  M,  the  order  of  the  fit.  This  is  the  case  because  the  increases  in 
actual  overhead  is  faster  for  the  batch  than  the  recursive  as  a  function  of 
the  window  length. 


4.4  Discussion  of  Results 


As  previously  shown,  the  recursive  method  provides  results  as  accurately 
as  the  batch  method  with  much  less  computational  speed  and  effort.  In  addi¬ 
tion,  the  recursive  method  provides  a  useful  estimate  even  before  the  data 
window  is  filled. 


The  theoretical  speed  advantage  shows  the  great  improvement  of  this 
method  over  the  batch  method.  The  actual  measured  speed  comparison  may  be 
made  to  move  toward  this  ideal  SR  by  optimizing  the  code  used  in  programming. 
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5.  SUMMARY  AND  CONCLUSIONS 

This  study  developed  an  efficient  recursive  algorithm  to  implement  a 
moving  average  Least  Mean  Square  Fit  (LMSF)  procedure.  The  following  briefly 
summarizes  the  significant  findings  of  this  study. 

1.  A  recursive  formulation  of  a  moving  average  LMSF  can  be  implemented 
with  a  theoretical  ratio  in  computation  speed  gain  (recursive  over 
batch)  of  (3M  +  4)/2N,  where  M  is  the  order  of  the  LMSF  and  N  is  the 
window  length. 

2.  It  was  shown  in  Section  3  that  recursive  formulation  is  identical 
to  the  batch  in  terms  of  filter  memory  length  and  data  storage 
requirement. 

3.  Ignoring  the  potential  singular  value  inversion  and  other  numerical 
problems,  recursive  moving  average  LMSF  gives  outputs  identical  to 
the  batch  LMSF  when  the  data  window  is  filled.  However,  prior  to 
attaining  the  full  data  window,  the  recursive  approach  has  the 
additional  advantage  that  it  could  also  provide  meaningful  outputs 
if  good  a  priori  information  is  used  to  initiate  the  recursion. 

4.  The  recursive  and  the  batch  LMSF  have  been  implemented  and  tested 
extensively  on  a  VAX  11/780  computer.  The  measured  performance  in 
terms  of  speed  and  output  behavior  agree  well  with  the  analysis. 
(Detailed  discussions  of  the  simulation  results  were  presented  in 
Section  4.) 

5.  The  significant  reduction  in  computational  load  for  the  moving 
average  LMSF  could  make  it  applicable  for  real-time  processing. 

6.  The  moving  average  LMSF  can  be  used  for  a  wide  variety  of  applica¬ 
tions  in  data  smoothing,  noise  filtering,  and  numerical  interpola¬ 
tion. 
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7.  Although  this  study  developed  the  moving  average  recursive  LMSF 
assuming  a  time  data  sequence,  the  result,  however,  is  applicable 
to  any  other  sequences  such  as  frequencies  and  wave  numbers. 

8.  We  have  shown  the  speed  advantage  of  the  recursive  approach.  How¬ 
ever,  we  have  ignored  the  relative  comparison  from  the  numerical 
point  of  view.  This  important  comparison  should  be  addressed  in 

a  future  study. 
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0001 

0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 

0011 

0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 

0021 

0022 

0023 

0024 

0025 

0026 

0027 

0028 

0029 

0030 

0031 

0032 

0033 

0034 

0035 

0036 

0037 

0038 

0039 

0040 

0041 

0042 

0043 

0044 

0045 

0046 

0047 

0048 

0049 

0050 

0051 

0052 

0053 

0054 

0055 

0056 

0057 


C  RECURSIVE  METHOD  FOR  FITTING  A  POLYNOMIAL  TO  DATA 

C  -UNIFORM  SAMPLING  ASSUMED- 

REAL*8  C(10,10) ,CM(8) ,YN(10) ,YO(10) , AN( 10 , 1 ) , Z ( 250) , 

ScH(  250 , 10 )  , HTET(  10,10)  , HTH(  10 ,  10 )  , HT(  10,250)  , HI (10, 10) 

IDGT=0 
WRITE  (6,3) 

3  FORMAT  (’  INPUT  SYSTEM  ORDER  (0<=  M  <=9)') 

READ  * , M 
WRITE  (6,6) 

6  FORMAT  ('  INPUT  THE  LENGTH  OF  THE  DATA  WINDOW  (l<  LW  <250)') 
READ  * , LW 

OPEN (UNIT=10,NAME=‘ OUT. DAT* ,TYPE='OLD' , FORM= ' FORMATTED ' ) 

DO  1=1, LW 
Z(I)=0 .0 
END  DO 

DO  1=1, M+l 
YO( I )=0 .0 
AN( 1,1 )=0 . 0 
END  DO 

C  CREATE  THE  H  MATRIX 

DO  1=1, M+l 
DO  K=1 , LW 

IF  (I-1.EQ.0)  THEN 
H(K,I)=1 
ELSE 

h(k, I ) =(K— 1 ) ** ( 1-1 ) 

END  IF 
END  DO 
END  DO 

C  CREATE  H  TRANSPOSE 

DO  1=1, LW 

DO  J=1 , M+l 

HT(J,I)=H(I,J) 

END  DO  " 

END  DO 

C  CREATE  ( HTH) ** ( —1 ) 

CALL  Ml MULT  ( HT , H, HTH, M+l , LW, M+l ) 

DO  1=1, M+l 
DO  J=1 , M+l 

HI ( I , J) =HTH( I , J) 

END  DO 
END  DO 

CALL  INVERT (HI, M+l) 

DO  1=1, M+l 
DO  J=1 , M+l 

HI ( I , J)=HI ( I , J) *LW 
END  DO 
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0058 

0059 

END  DO 

0060 

0061 

C 

CALCULATE  MATRIX  OF  BINOMIAL  COEFFICIENTS 

0062 

IF  (M.GT.l)  THEN 

0063 

DO  1*2, M 

0064 

C(l-l,l)=l 

0065 

END  DO 

0066 

DO  L=2 , M 

0067 

C(L— 1 , L+l ) *1 

0068 

Kl=l 

0069 

DO  1*2, L 

0070 

K1=K1*I  1  Kl= ( M— 1 )  1 

0071 

END  DO 

0072 

DO  J=2,L 

0073 

K2=l 

0074 

K3=l 

0075 

DO  N=1 ,  ( L- J+l ) 

0076 

K2=K2*N  1  K2= ( M-I ) 1 

0077 

END  DO 

0078 

DO  N=2 , J 

0079 

K3=K3*(N-1)  1  K3=( 1-1)1 

0080 

END  DO 

0081 

C(L-l,J)=Kl/K2/K3 

0082 

END  DO 

0083 

END  DO 

0084 

END  IF 

0085 

0086 

0087 

C 

CALCULATE  CM  VECTOR 

0088 

IF  (M.GT.l)  THEN 

0089 

DO  1=3, M+l 

0090 

DO  J«1,I 

0091 

CM{ 1-2 )=CM( 1-2 )+( LW-1 ) ** ( J-l ) *C( 1-2 , J) 

0092 

END  DO 

0093 

CM( I— 2 )=CM( 1-2 ) / FLOAT ( LW) 

0094 

END  DO 

0095 

END  IF 

0096 

0097 

0098 

C 

BEGIN  POLYNOMIAL  MATCHING  RECURSION 

0099 

IP=1 

0100 

0101 

L=0 

0102 

DO  WHILE(L.LT.l) 

0103 

0104 

0105 

C 

READ  IN  NEW  DATA 

0106 

READ  (10,*)  ZA 

0107 

0108 

WRITE  (9,*)  EST , ZA 

0109 

0110 

IF  (ZA.EQ.0.)  GOTO  100 

0111 

YN ( 1 ) =YO( 1 )+( ZA-Z ( IP )) /FLOAT ( LW) 

0112 

IF  (M.GT.0)  THEN 

0113 

YN( 2 )=YO( 2 )+ZA-YN( 1 ) 

0114 

IF  (M.GT.l)  THEN 
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0115 

DO  1=3, M+l 

0116 

R=0 

0117 

DO  J=1 , 1-1 

0118 

R=R+C(l-2,J)*YN(J) 

0119 

END  DO 

0120 

YN ( I )=YO( I )— R+CM( I— 2 ) * 

0121 

END  DO 

0122 

END  IF 

0123 

END  IF 

0124 

0125 

C 

CALCULATE  POLYNOMIAL  COEFFICIENTS 

0126 

0127 

CALL  M3 MULT  (HI , YN, AN, M+l , M+l , 1 

0128 

0129 

C 

CALCULATE  ESTIMATE  OF  NEXT  POINT 

0130 

0131 

EST=0 

0132 

DO  J=1 , M+l 

0133 

EST=EST+AN( J , 1 ) * ( LW+1 ) ** (J-l 

0134 

END  DO 

0135 

0136 

C 

UPDATE  DATA  WINDOW 

0137 

0138 

Z(IP)=ZA 

0139 

IP=IP+1 

0140 

IF  (IP.GT.LW)  IP=1 

0141 

0142 

C 

MAKE  CURRENT  RECURSION  RESULTS  'OLD 

0143 

0144 

DO  1=1, M+l 

0145 

Y0(I)=YN(I) 

0146 

END  DO 

0147 

END  DO 

0148 

0149 

100 

STOP 

0150 

END 
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